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We argue that Josephson junction networks may be engineered to allow for the emergence of 
new and robust quantum coherent states. We provide a rather intuitive argument showing how the 
change in topology may affect the quantum properties of a bosonic particle hopping on a network. 
As a paradigmatic example, we analyze in detail the quantum and thermodynamic properties of 
non-interacting bosons hopping on a comb graph. We show how to explicitly compute the inhomo- 
geneities in the distribution of bosons along the comb's fingers, evidencing the effects of the topology 
induced spatial Bose-Einstein condensation characteristic of the system. We propose an experiment 
enabling to detect the spatial Bose-Einstein condensation for Josephson networks built on comb 
graphs. 



I. INTRODUCTION 



Quite recently it has been evidenced that quantum states may support new kind of orders, which cannot be 
characterized by broken symmetries and, thus, cannot be described by the conventional Ginzburg-Landau theory 
0. Quantum orders may be viewed as the pertinent description of the pattern of the quantum entanglement in a 
quantum many-body ground-state. Quantum or topological orders have already been extensively used in the analysis of 
Fractional Quantum Hall systems 0, leading to an elegant explanation of their robust (against a weak but, otherwise, 
arbitrary perturbation) ground-state degeneracy and evidencing the intimate connection between the ground-state 
and the statistics of quasi-particles. The robustness of the ground-state degeneracy makes quantum ordered states 
very relevant candidates for the engineering of quantum devices naturally taming the intrinsic decoherence of other 
solid state devices 0). Very recent studies hint to the new and exciting possibility that the topology of Josephson 
junction networks (JJN) may be crucial for inducing novel and unexpected macroscopic quantum phenomena y| as 
well as for opening |5j to the possibility for an explicit realization of a quantum ordered state. 

In this paper we shall focus our attention on the recently discovered possibility that JJN built on suitable graphs may 
support new and interesting quantum macroscopic states induced only by a pertinent engineering of the geometry and 
topology of the graph supporting the network HQ. In particular, we shall consider JJN built on comb-like graphs, 
evidencing the fact that - already for this very simple graph topology - the array exhibits quantum macroscopic 
coherence at low temperatures. After discussing the theoretical possibility of a spatial Bose-Einstein condensation 
(BEC) of Cooper pairs for a JJN built on a comb-graph, we shall propose here a simple experiment which could 
enable to observe BEC in these systems. 

A network of superconducting classical Josephson junctions is usually described by the Hamiltonian 0: 



Hjj = -J ^2 cos( 

(x,y;x',p') 



'x,y 



(1) 



On each site of the network there is a superconducting grain and the junctions are located between neighboring sites 
with Josephson energy Jq; (j) XiV is the phase of the superconducting order parameter at site labeled by (x,y). The 
parameter Jo is connected to the critical Josephson current I c by the relation I c = 2eJo/h [8j. The sum (x,y; x',y') 
runs over all the distinct nearest-neighbors pairs. The Hamiltonian Q is closely related to the Hamiltonian describing 
non-interacting bosons hopping on a generic network 

H = -t Yl (4yV,y' + (2) 

(x,y; x',y') 

In Eq. @, t is the hopping parameter while a x y (a XyV ) is the creation (annihilation) operator for bosons; h x ,y — 
a x y a Xj y is the number operator at site (x, y). The filling, i.e., the average number of particles per site, is defined as 
/ = Nt/L 2 , where Nt is the total number of bosons and L 2 is the number of sites. When / S> 1 and the fluctuations 
of the particle numbers per site are much smaller than /, one can safely substitute the operator a x ,y with */ ' N x ^ y e l ^ :c ' y , 
where N XiV is the number of particles at the site (x,y) 0, El- As a result, the Hamiltonians and may be 
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FIG. 1: Comb graph. 



regarded as equivalent, provided that 

Jo « 2ty/N XiV N x ,, y , « 2t/. (3) 

The plan of the paper is the following: in Sec. II we analyze the spectrum of bosons hopping on a comb-shaped 
graph. There we shall also use a variational approach aimed to elucidate how a mere change in the network's topology 
is crucial for explaining the emergence of localized states in the single-particle spectrum of free bosons hopping on 
a pertinent graph. In Sec. Ill we analyze the thermodynamic properties of non-interacting bosons hopping on a 
comb-graph: we compute the critical temperature T c and the inhomogeneous spatial distribution of bosons over the 
array as a function of the scaled temperature T/T c . In Sec. IV we propose a method to detect BEC in a comb-shaped 
JJN: for this purpose we shall determine the Josephson critical current I c as a function of the scaled temperature 
T/T c and of the location of the junction. 

II. TOPOLOGY INDUCED LOCALIZATION FOR BOSONS HOPPING ON A GRAPH 

In this section, we shall illustrate how a mere change in the topology of a lattice may affect the spectrum of 
a quantum mechanical system by analyzing the properties of non-interacting bosons hopping on a comb-shaped 
network. 

A comb graph (see Fig. ^) is composed of one-dimensional chains (fingers) grafted periodically on a linear chain 
(backbone). Each site of the comb can be naturally labeled introducing two integer indices [x,y), where x labels the 
different fingers and y represents the distance from the backbone. Each site on the finger is linked to two neighbors 
whereas each site of the backbone has four neighbors. The Hamiltonian (J2J) on a generic graph can be written as: 

H = —t y ' A x ,y; x',y' a x,yG>x',y' ■ (4) 
x,y ; x' ,y' 

The topology of the network is fully described by the adjacency matrix A x ,y, x ',y' which is equal to 1 if (x, y; x', y') is 
a link and otherwise. The single-particle energy spectrum a on a L x L comb array may be found by solving the 
eigenvalue equation 0,0: 

-t ^2 A x>v . x ,y ip E (x', y') = Eip E (x, y). (5) 

x',y' 

The spectrum a is formed by L 2 states and it is divided in four regions: {E }, <j_, o~q and a + jj. E is the ground-state 
energy which, with periodic boundary conditions and in the thermodynamic limit L — > oo, is given by E = -2y/2t. 
<7o is the part of the spectrum corresponding to delocalized states with energies between — It and 2t; the density of 
states of do is 



P[E) = V7W=W ' 



(6) 



FIG. 2: From square lattice to comb. 



just as for a particle hopping in a linear chain. The so-called hidden spectrum |4| is given by the union of cr_ = 
{— 2tyl + cos 2 (k)} and a + = {2t^/l + cos 2 (k + n/2)}, where k = 2im/L is the wave vector along the backbone 
introduced by the Fourier transform of the x coordinate and n = 1, . . . , (L — l)/4; the hidden spectrum does not 
contribute to the normalized density of states in the thermodynamic limit since the number of these states is L 4] . The 
ground-state and the eigenstates belonging to <j_ and a + - due to the particular topology of the array - are localized 
along the backbone and exhibit an exponential decay in the direction of the fingers. Omitting the normalization 
factor, the ground-state wavefunction is given by 3 

V^ mb (z,y) = e- M/ « , ( 7 ) 

where £o satisfies the transcendental equation sinh(l/£o) — coth[L/(2£o)] = 0; the associated eigenvalue is Eq = 
— i(e~ 1 /'> + e 1 '^ ). As we shall see in the next section, bosons are allowed to spatially condense in this localized state. 



It is instructive to provide a more intuitive path to evidence how a mere change in the topology of a network may 
give rise to the localized ground-state for bosons hopping on a comb graph. For this purpose one may start from a 
square lattice and imagine to systematically cut the bonds linking the sites of the square lattice in the x-direction at 
a chosen value of the y coordinate (see Fig. if one cuts all the bonds in the x direction but the ones at the origin 
y = 0, one obtains the comb lattice. Starting from a square lattice with (2M + 1) x (2M + 1) sites, the procedure 
to follow is to remove all the links in the x direction at distance \d\ = M from the origin and then gradually cut 
stripes of links at \d\ = M — 1, . . . , 1. As a result, one gets M different arrays, each described by a different adjacency 
matrix A™ yl x , y , , where 2m is the number of cuts and [2M + 1) — 2m = 2b + 1 is the number of remaining backbones 
(6 = M — m). For the square lattice one has: 

A S x q y - x >^> = {5x,x' + l + Sx,x'-l) ■ 5y,y> + 5x,x> ' (*y, a ' + l + <W'-l)> ( 8 ) 

for a square with 2m cuts the adjacency matrix is, for m > 1, given by 

m— 1 

A 7,y; x>,y> = A 7,y;x> , y> ~ 2 (<W+1 + S x,x'-l)Sy,y'6y tM - q - (9) 

q=0 

To determine the energy of the ground-state as a function of the number of cuts, it is most convenient to use a 
variational approach. For this purpose one may use a trial ground-state wavefunction ip^(x,y) and then minimize 
the value of the energy E(£,b;M) with respect to the parameter £. With periodic boundary conditions, ip^(x,y) 
should interpolate - as m increases - between the ground-state wavefunction corresponding to bosons hopping on a 
square lattice (a constant with associated eigenvalue Eq = —At) and the wavefunction J7J, which is the ground-state 
wavefunction for bosons hopping on a comb graph. 

The simplest choice of the trial ground-state wavefunction taking into account both limits is given by 



iP 6 (x,y) = e-\y\/t, 



(10) 
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FIG. 3: Variational energy 1131k in units of t, of the ground-state of the single-particle eigenvalue equation © as a function of 
the number of backbones b in the thermodynamic limit M — > oo. When 6 = 0, one finds the ground-state energy of the comb 
graph E = -2\/2t. 



where £ > is the variational parameter. The variational energy of the ground-state of bosons hopping on a lattice 
with 2m stripes of links removed is given by: 

M M 

E^m;M) = t = t — t^L ^ . 

mm) H, . , . , 

x,y=—M 

As a function of the number of backbones, one has: 



e 



-2b n 



| 4e (i+26)/€ + e 2M/i \2 _ e -26/4 ( X + e i/?) 2 ] I 



^.6^0=* e 2M/c + e 2(M + i)/; _ 2 ( 12 ) 

If one minimizes Eq. I|12[l with respect to the parameter £, one obtains the variational ground-state wavefunction and 
its associated energy as a function of b. 

The results are summarized in Figs. EE and are compared with the results of a numerical evaluation of the ground- 
state wavefunction and its energy. In Fig. [3]we plot the variational energy in the thermodynamic limit which is taken 
by keeping b (i.e., the number of remaining backbones) constant and performing the limit M — ► oo; the energy is 
given by 
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2e - 2b /« - (l + e 1 /?) 



£(6,0=*— r—7T? -■ (!3) 



2 2/« 

As expected, effects due to the network topology clearly emerge if the number of remaining backbones is a zero- 
measure set with respect to the number of stripes in a square lattice. In Fig. 0|we plot the numerical ground-state 
wavefunction of bosons hopping on a network with (51 x 51) sites for different values of m. Figure 0] evidences that the 
simplest trial wavefunction given in Eq. (|l(Jf> is rather accurate in describing the shape of the exact wavefunction for 
a number of cuts m ~ M . In Fig. [S]we compare the energies pertinent to the ground-states discussed in this section: 
crosses correspond to the numerical results and the solid line to the ansatz <|1UI) . We also plot, as a function of the 
number of cuts m, the energy of the ground-state of bosons hopping on a square lattice: for m — the energy is — 4t, 
while for m = M (corresponding to the comb lattice) is — 2t. The dashed line is the energy of the ground-state Q of 
bosons hopping on a comb lattice as a function of the number of cuts m: for m — M the energy is — 2\/2t, while for 
m = the wavefunction (JJJ is, of course, an excited state. Figure |S] clearly illustrates how the localized ground-state 
of bosons hopping on the comb becomes energetically more favorable when one increases the number of cuts. 

The variational and numerical analysis carried out above provides us with an intuitive explanation of the appearance 
of localized states in the spectrum of a quantum particle hopping on a pertinent graph; furthermore, it clearly relates 
this property to the topology change of the network. 
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FIG. 4: The normalized single-particle ground-state wavefunction i/)(x,y) computed numerically as a function of the distance 
y from the backbone for a comb graph with 51 x 51 sites and different numbers of cuts m = 25, 20, 15, 10, 5, and 0. m — 
corresponds to the square lattice and m = 25 to the comb. 
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FIG. 5: Single-particle ground-state energies(in units of t) as a function of the number of cuts m for a 51 x 51 lattice: crosses 
correspond to the numerical results, the solid line to the ansatz < I IH . We also plot the energies of the ground-state of bosons 
hopping on a square lattice (dotted line) and on a comb lattice (dashed line). 



III. THERMODYNAMICS OF BOSONS HOPPING ON A COMB GRAPH 



The thermodynamic properties of non-interacting bosons hopping on a comb graph hint to the possibility of a 
topology induced spatial BEC even if d < 2 0, To elucidate this new phenomenon, it is most convenient to 
introduce the macrocanonical ensemble to determine the fugacity z = ef-^ E ^ as a function of the temperature 
for a lattice with L x L sites, the equation determining z is given by 

fL> = N T = r ■ -gg L : ■ (14) 



EEcr 



z -i e 0(E+V5t) 



In Eq. H14|). d(E) is the degeneracy of each single-particle eigenstate of the Hamiltonian Q and (3 = l /fepT; the sum 
in Eq. I|14|) is over the entire spectrum a. For free bosons hopping on a comb graph, one has 0, 0, Il2j | 

L 2 P (E) 

EEa Z-'-et 

where Ne (L,T), N a _(L,T) and N a+ (L,T) denote, respectively, the number of particles at a certain temperature 
T in the ground-state and in the two regions er_ and cr+ of the hidden spectrum. p(E), with E G do, is the energy 



N T = N Eo (L,T)+N a _(L,T)+N a+ (L,T)+ I dE — ^J t) _ , (15) 
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density of states defined in Eq. ©, with Co being the region of the spectrum corresponding to delocalized states. It 
is also useful to define the number of particles per site in each part of the s pec trum as ue = Ne / L 2 , n aa = N ao /L , , 
Tifj = N a _/L 2 and n a+ — N a+ /L 2 . In the thermodynamic limit, one has [12] 

1 1 



n Eo (T) = lim 



L^oo L 2 z- 1 -!' 



2 {L ~ 1]/i 1 
n a _ (T) — lim - \ 



L 2 7 _l p /3i[v / 8-2 x /l+cos 2 (27rn/L)] 
n— 1 * e v 

and 

(£-l)/4 



L 2 ^ r -l„/3i[V8+2 x /l+cos 2 (27rn/L+7r/2)] _ -i 
n— 1 ^ c v J- 

< lim I- -V- =0 V T. 

Thus, in the thermodynamic limit, a+ is not macroscopically occupied at any temperature and does not play any role 
in describing the thermodynamics of the system. 

The last term of the right-hand side of Eq. (|15|l is the number of bosons in the delocalized (chain-like) states. 
The presence of the hidden spectrum changes the behavior of the integral evaluated in the interval {—It, 2t}, since it 
reduces it to the one describing non-interacting bosons on a linear chain with an impurity in one of the sites. As a 
result, letting z — > 1, the integral converges even at finite temperatures making possible a spatial BEC in d < 2. 



A. Critical temperature and condensate fraction 

If one defines T c as the critical temperature at which BEC occurs, for any T < T c the ground-state is macroscopically 
filled. Since, at the critical temperature, riE (T c ) — n a± {T c ) = 0, the equation allowing to determine T c as a function 
of the parameters / and t reads: 

r , dE * . as) 

J -it \/4t 2 - E 2 e {E+^t)/(k B T c ) _ i v ; 

Equation (|16f) can be solved numerically for any value of /. When / 3> 1, one may expand the exponential in Eq. 
(|16|l to the first order in the inverse of the critical temperature T c getting 

fnTdE - 1 ^ (17) 

J-2t 7T V4t 2 - E 2 E + V&t 2t 

The critical temperature T c is then a linear function of both t and / given by: 

T c ^. (18) 

Equation Ijl8(l has been checked numerically and it is in excellent agreement with the numerical solution of Eq. Ijlfijl 
for / ^> 1, the error being of order 1//. By means of Eq. ©, one finds that the critical temperature for the occurrence 
of BEC is 

T c «^. (19) 
kb 

One may now use Eq. (|18|) to determine the condensate fraction as a function of the scaled temperature T/T c . In the 
thermodynamic limit, the number of particles in the delocalized states is given by 

»- fe) " SS. L " L, "' £ »^oTT » *£■ <»> 
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In the last equation the exponential has been expanded to the first order in [3: this approximation holds for / 3> 1 
and it is in very good agreement with the numerical evaluation of the integral in Eq. I(2U() also in a large neighborhood 
below T c . From Eqs. (fT5)l and lf277jl one gets the number of particles in the localized states N Q = Ne + N a _ : the 
fraction of condensate, for T <T C , is then given by 

For / ranging from 10 3 to 10 9 , the results provided by Eq. (|21|l differ from those obtained by the numerical evaluation 
of Nq from Eq. 1)15(1 by less than 1%. Equation ((21(1 clearly shows that the condensate has dimension 1; cig ar-shaped 
one-dimensional atomic Bose condensates support, in fact, a condensate fraction given by Eq. (j2H> [Hill- 



B. Distribution of bosons along the fingers 



In the following we shall determine the distribution of the bosons over the comb graph, since - due to the spatial 
condensation - one expects an inhomogeneous distribution of the bosons along the fingers. The average number of 
bosons Nb{x, y) on a site (x, y) does not depend on x - due to the translational invariance along the backbone - but 
only on the distance y from the backbone, and - at any temperature T - is given by: 



N B (y;T/T c ) = N Eo (T/T c ) |Ve (j/)| 2 + ]T N a _ (E n ;T/T c ) \^ E Jy)\' 

1 



L 



dE p{E) 



,P(E+\/8t) 



\Mv)\' 



(22) 



In Eq. ((22(1 ipE (y) is the wavefunction corresponding to the ground-state of the single-particle spectrum and ipE n (y) 
are the eigenfunctions corresponding to the energies E n of the hidden spectrum N a _ (E n ) is the number of particles 
with energies E n € cr_ and Ne„ is the number of particles in the ground-state. In the last term of Eq. 1(22(1 ipE(y) 
are the delocalized eigenfunctions of the eigenvalue equation J5|. 

To determine Ng (y;T/T c ), one needs an expression for Ne and N a _. For this purpose, it is useful to define the 
scaled temperature 



T 

T= T- 

For r < 1 and in the thermodynamic limit, N a _ (E n , r) is given by ^ 

N a _{E n -r) = lim L 2 



L — ►oo 



Ne (t) 



and it depends on the number of particles in the ground-state Ne . Using the fact that ksT c 
follows that the number of particles in the hidden spectrum is given by 



N„_ (r) = J2 N *- ( E n) = ~ N Eo + N Eo 1 



Nn 



V2N Eo 



coth 



Nn 



V2N Eo 



(23) 
(24) 

2tf, from Eq. ^ it 
(25) 



Requiring that N = N a _ 
the scaled temperature r: 



Ne and using Eq. ((21(1 . one finds the way to determine Ne /Nt as a function only of 



JVt(I-t) = N Eo < 



N; 



y/2N Eo 



coth 



N, 



(26) 



V2N Eo 

allows for an exact numerical evaluation 



Solving Eq. I(2tj|) and substituting back the value obtained for Ne in Eq 
ofEq. 

Conventional wisdom supported by numerical evidence suggests however that - apart from a small range of tem- 
peratures near T c - the largest contribution to Nq comes from Ne - Thus, it is physically appealing to assume the 
following form for Ne , namely, 



Neo = N T (1 - t) 9 (t). 



(27) 
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The exponential behavior typical of the eigenfunctions corresponding to the localized states leads, for r < 1, to an 
increase of Nb{d',t) on the sites near the backbone while, when y 3> 1, the behavior is dominated by the last term 
in the right-hand side of Eq. Ij32(l . Thus, away from the backbone, once the filling is fixed, N B depends only on the 
scaled temperature r and it is given by 

N B (y;T/T c ) _ T 

J r = -. (33) 

In the next section we shall show how this property may be useful to compute the observable effects induced by the 
existence of a spatial BEC in condensed matter systems. 



IV. FURTHER OUTLOOK: SIGNATURE OF BEC FOR JOSEPHSON NETWORKS ON A COMB 

We proposed so far a theoretical framework establishing the role played by the network's topology in understanding 
the quantum and thermodynamic inhomogeneities arising in a system of bosons hopping on a comb graph; the 
approach enables one to obtain a simple and explicit analytical expression for the inhomogeneous distribution of 
bosons along the comb's fingers signaling the dramatic effect of the topology induced spatial BEC evidenced in Refs. 
] . With little modifications our analysis could be carried out also for diverse network's topologies supporting BEC 



We shall now evidence how to detect the signature of BEC for a JJN built on a comb graph. For this purpose, 
one may think to perform a measurement of the I-V characteristic of a single finger of the JJN and of the critical 
Josephson current along the finger. If one feeds, in fact, an external current I ex t at the extremities of the finger, one 
should expect to observe no voltage unless I ex t is larger than the smallest of the critical currents of the junctions 
along the finger. Since, below T c , the critical Josephson current of the finger is given by the smallest of the critical 
currents of the junctions positioned along the fingers, the measurement of the I-V characteristic of the finger should 
provide a measurement of the critical current of the junctions at the top of the fingers (i.e., at y 3> 1). 

To make a definite prediction, one needs to estimate the value of the Josephson critical current of a single junction 
of the network as a function of both the temperature and the distance from the backbone. Above T c the Cooper pairs, 
for a large array (L 3> 1), are uniformly distributed over the network: the Josephson critical current is the same at 
each junction 

It (y) « t (Ae/h) ■ VN B (y+l)N B (y) « J • (2e/ft) (34) 

where 7^ is the Josephson critical current above T c . According to Eq. the Josephson critical current depends 
only on the position and on the population of the sites. 

The relation between the Josephson critical current of the junction between the sites (x, y) and (x, y+ 1) above and 
below the critical temperature T c does not depend on x and it is given by 



Ic(y,r) _ y/N B (y + l,T)N B (y,T) ( 

~im i ' (35) 

where lf(y;r) is the Josephson critical current below T c at distance y from the backbone. N B (y,r) is given by Eq. 

(|32[) . The result l|35l) - for a given / - is independent from the total number of bosons in the system. Far from the 
backbone (y 3> 1), one simply has 

Ic(v,r) „ N B (y, T ) 

~^?W 7 ' (36) 

and, using Eq. one gets 

If (V5 r) T 



i c A (y) t c 



(37) 



Ic(y;r) at the top of the finger can be experimentally measured as the critical Josephson current of the finger. In 
Fig □ we plot Eq. j37|). 

BEC in a JJN built on a comb graph predicts then a rather sharp decrease of the Josephson critical current for a 
junction located away from the backbone and this behavior affects the measurement of the I-V characteristic along a 
given finger of the JJN. The linear dependence exhibited by the solid line in Fig. [7|is consistent with the observation 
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FIG. 7: Josephson critical current as a function of T/T c computed at distances from the backbone y ^> 1. if (y;T/T c ) is in 
units of lf(y) (the critical current of the junctions above T c ) and is therefore equal to 1 for T > T c . 



that, in this system, the condensate has dimension 1. We stress that this phenomenon should be already observed for 
a classical JJN described by the Hamiltonian Q. 

The slope of the linear plot in Fig. [7|provides a direct estimate for the fraction of condensate N /Nt ~ 1 — ksT/J 
for a JJN built on a comb graph. In fact, far away from the backbone (y 3> 1), from Eqs. I)2ip. I|3(jfl and (|37[) . one 
has that the fraction of condensate is given by 1 — \lf(y; t)/ if . 

To conclude this section, we notice that a similar analysis applies to a bosonic gas in comb-shaped deep optical 
lattice, provided that one defines the Josephson current of a single bosonic Josephson junction between two neighbor 
sites at y and y + 1 as I c w 2t^/ 'TV 'g(y)N s{y + 1) (where Nb(v) is the number of bosons in y). Above T c , all the 
Nb(u) will be equal to the filling / at the equilibrium and if = 2tf. Below T c , the equilibrium values for Nb(v) 
change along the fingers according to Eq. (1221) and the ratio between lf(y,r) and if gives, for y 3> 1, directly the 
fraction of noncondensate atoms. 
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